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Abstract. 

We investigate a microscopic motor based on an externally controlled two-level 
system. One cycle of the motor operation consists of two strokes. Within each stroke, 
the two-level system is in contact with a given thermal bath and its energy levels are 
driven with a constant rate. The time evolution of the occupation probabilities of the 
two states are controlled by one rate equation and represent the system's response 
with respect to the external driving. We give the exact solution of the rate equation 
for the limit cycle and discuss the emerging thermodynamics: the work done on the 
environment, the heat exchanged with the baths, the entropy production, the motor's 
efficiency, and the power output. Furthermore we introduce an augmented stochastic 
process which reflects, at a given time, both the occupation probabilities for the 
two states and the time spent in the individual states during the previous evolution. 
The exact calculation of the evolution operator for the augmented process allows us 
to discuss in detail the probability density for the performed work during the limit 
cycle. In the strongly irreversible regime, the density exhibits important qualitative 
differences with respect to the more common Gaussian shape in the regime of weak 
irreversibility. 
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1. Introduction 

Non-equilibrium phenomena in the presence of time-varying external fields are of 
vital interest in many areas of current research [H l2l 13]. Examples are aging and 
rejuvenation effects in the rheology of soft-matter systems and in the dynamics of 
spin glasses, relaxation and transport processes in biological systems such as molecular 
motors, ion diffusion through membranes, or stretching of DNA molecules, driven 
diffusion systems with time-dependent bias, and nano-engines. With minimization of 
the system size thermal fluctuations become increasingly relevant. In these systems 
it is useful to introduce microscopic heat and work quantities as random variables 
whose averages yield the common thermodynamic quantities. Averages over functions 
of these microscopic heat and work quantities yield generalized fluctuation theorems 
[H El El [71 El El [ini [m [El [13 [E] . in this context mesoscopic engines operating between 
different heat baths under non-equilibrium conditions have received increasing attention. 
The variety of models can be roughly classifled according to the dynamical laws involved. 
In the case of the classical stochastic heat engines, the state space can either be discrete 
or continuous (c.f., for example, [151 [13 [12] and the references therein). Examples of 
the quantum heat engines are studied, e.g., in [T8[ 119] . 

The traditional consideration of efficiency of heat engines operating between two 
baths at temperatures Ti and T2 leads to the Carnot upper bound rjc = 1 — T1/T2. 
The bound is only achieved under reversible conditions where the state changes require 
infinite time and hence the power output is zero. Real heat engines generate a finite 
power output Pout = W^out/^cycie, i-G., they perform work PVout during a cycle of a finite 
duration tcycie- Thus an appropriate way to characterize the engines is to compare 
their efficiencies at maximum power. On the macroscopic level this quantity is roughly 
bounded by the Curzon-Ahlborn value rjcA = 1 — V^V^ PO]. Alternative expressions 
for quantifying efficiency have been discussed [17] which are based on mean quantities, 
e.g., on the mean work done during the operational cycle. On the mesoscopic level, the 
work is inherently a fluctuating quantity and one should be able to calculate not only 
its mean value but also its fluctuation properties. 

In this paper we study a simple model of mesoscopic heat engines operating between 
two different heat baths under non-equilibrium conditions. The working medium consist 
of a two- level system. The cycle of operation includes just two isothermal branches, 
or strokes. Within each stroke, the system is driven by changing the energies of the 
two states and we assume a constant driving rate, i.e. a linear time-dependence of the 
energies. The response of the working medium is governed by a master equation with 
time-dependent transition rates. The speciflc form of the rates guarantees that, provided 
the two energies were flxed, the system would relax towards the Gibbs equilibrium state. 
Of course, during the motor operation, the Gibbs equilibrium is never achieved because 
the energies are cyclically modulated. At a given instant, the system's dynamics just 
reflects the instantaneous position of the energy levels. After a transient regime, the 
engine dynamics approaches limit cycle with the periodicity of the driving force. We will 
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focus on the properties of this hmit cycle. In particular, we calculate the distribution 
of the work during the limit cycle. 

Our two-isotherm setting imposes one important feature which is worth 
emphasizing. As stated above, at the end of each branch we remove the present bath 
and we allow the thermal interaction with another reservoir. This exchange of reservoirs 
necessarily implies a finite difference between the new reservoir temperature and the 
actual system (effective) temperature. Even if the driving period tends to infinity, we 
shall observe a positive entropy production originating from the relaxation processes 
initiated by the abrupt change of the contact temperature. Differently speaking, our 
engine operates in an inherently irreversible way and there exists no reversible limit. 

The paper is organized as follows. In section [2] we solve the dynamical equation for 
the externally driven working medium. For the sake of clarity we first give the solution 
just for an unrestricted linear driving protocol using a generic driving rate and a generic 
reservoir (section l2.ll) . Thereupon, in section [2^ we particularize the generic solution 
to individual branches and, using the Chapman-Kolmogorov condition, we derive the 
solution for the limit cycle. In section [3] we employ the recently derived pi] analytical 
result for the work probability density under linear driving. Again, we first give the 
result for the generic linear driving and then we combine two such particular solutions 
into the final work distribution valid for the limit cycle. The results from section [2] and 
section [3] enable a detailed calculation of the energy and entropy flows during the limit 
cycle in section H] and allow for a discussion of the engine performance in section [51 

2. Description of the engine and its limit cycle 

Consider a two-level system with time- dependent energies Ei{t), i = 1,2, in contact with 
a single thermal reservoir at temperature T. In general, the heat reservoir temperature 
T may also be time-dependent. The time evolution of the occupation probabilities Pi{t), 
i = 1,2, is governed by the master equation [22] with time-dependent transition rates 
specified by the reservoir temperature and by the external parameters. To be specific 
the dynamics of the system is described by the time inhomogeneous Markov process 
D{t) assuming the value i, i = 1,2, if the system resides at time t in the ith state. 
Explicitly, the master equation reads 

where I is the unity matrix and M(t, t') the transition matrix with elements Rij{t \ t') = 
{ i I R(t 1 1') i,j = 1, 2. These elements are the conditional probabilities 

R,,{t 1 1') = Prob { D(t) = 1 1 D(t') = J } . (2) 

If we denote by 0(t') the initial state at time t' with the occupation probabilities 
Pi{t') = {i \(j){t')), the occupation probabilities at the observation time t are described 
by the column vector \p{t,t')) = R(t 1 1') |0(t'))- 
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Due to the conservation of the total probability the system ([T]) can be reduced to 
just one non-homogeneous linear differential equation of the first order. Therefore the 
master equation ([T]) is exactly solvable for arbitrary functions \i(t), X2(t). The rates 
are typically a combination of an attempt frequency to exchange the state multiplied 
by an acceptance probability. We shall adopt the Glauber form 

AlW = T- r or/,r,f,, r,f,,n ^ = Ai (t) eXp { [i^i (t) - W] } , (3) 

1 + exp {-f3{t) [Ei{t) - E2{t)\} 

where z/"^ sets the elementary time scale, and /3(t) = l/kBT{t). The rates in equation ([1]) 
satisfy the (time local) detailed balance condition. 

The general solution of the master equation ([T]) for the transfer rates reads 

m\t')=I-l(_\ "M{l-exp[-Kt-t')]} + ^( "J ~\ U(t,0, (4) 



(5) 



where 

at,t') = 1^1^' dr exp [-u{t - t)] tanh [E,{t) - E^ir)]^ . 

The resulting propagator satisfies the Chapman- Kolmogorov condition 

R{t 1 1') = R{t I t")R{t" 1 1') (6) 

for any intermediate time t". Its validity can be easily checked by direct matrix 
multiplication. The condition simply states that the initial state for the evolution in the 
time interval [t",t\ can be taken as the final state reached in the interval [t',t"]. This 
is true even if the parameters of the process in the second interval differ from those 
in the first one. Of course, if this is the case, we should use an appropriate notation 
which distinguishes the two corresponding propagators. This procedure will be actually 
implemented in the paper. Keeping in mind this possibility, we shall first analyze the 
propagator for a generic linear driving protocol. 

2.1. Generic case - linear driving protocol 

Let us consider the linear driving protocol Ei{t) = h + v{t — t'), and E2{t) = —Ei{t), 
where h = Ei{0) denotes the energy of the first level at the initial time t', and v is the 
driving velocity (energy change per time). The rates ([3]) can then be written in the form 

^'W-^TT ^orT^- A.(f)^. °°-Pl-"'' -''>' . (7) 

1 + cexp[— — t )\ 1 + cexp[—\l{t — t')\ 

where Q = 2(3\v\ is the temperature-reduced driving velocity, and c = exp{—2f3h\v\/v) 
incorporates the initial values of the energies. 

Under this linear driving protocol one can evaluate the definite integral in (jlj) 
explicitly and rewrite the propagator as 

M(t|t')=I- ( _J Q ) {1-exp [-Kt-t')]}+ ( _J _[)l{t,t'), (8) 
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where 



, , , , exp (-fir) 
vc I dr exp [— z/(t — r)J 



nt 



1 + cexp (—fir) 
exp [(a — l)r] 



L/"l+ cexp (-r 



(9) 



Here we have introduced the dimensionless ratio a = u/Q of the attempt frequency 
characterizing the time scale of the system's dynamics and the time scale of the external 
driving, respectively. Naturally, this ratio will describe the degree of irreversibility of 
the process. Depending on the value a G (0, oo), the explicit form of the function •yit, t') 
reads 

' ^ exp {-ut) {exp [(a - l)fit'] 2Fi(l, 1 - a; 2 - a; -cexp (-fit')) 
-exp [(a - l)fit] 2Fi(l,l -a;2-a;-cexp(-fit))} 

1 + cexp (—fit) 



cexp (—fit) 



fi(t - t') 



In 



aG (0,1), 
a = l, (10) 



1 + cexp (-fit') J 
2Fi(l,a;l + a;-lexp (fit)) 

- exp [-afi(t - t')] 2Fi(l, a; 1 + a; -\ exp (fit')) , a > 1 , 

where 2Fi(a,/?;7; ■) denotes the Gauss hypergeometric function [23] . 

2.2. Piecewise linear periodic driving 

We now introduce the setup for the operational cycle of the engine under periodic 
driving. Within a given period, two branches with linear time-dependence of the state 
energies are considered with different velocities. Starting from the value hi, the energy 
Ei{t) linearly increases in the first branch until it attains the value /i2 > h, at time t+ 
and in the second branch, the energy Ei{t) linearly decreases towards its original value 
hi in a time t_ (see figure 1). We always assume E2(t) = —Ei{t), i.e. 

r hi + l^^^t, te[0,u], 

Elit) = -E,{1) = { (11) 

[ /i2 - ^^^^ (t - 1+) , te[u,u + t.] . 

This pattern will be periodically repeated, the period being tp = t+ + 1_. 

As the second ingredient, we need to specify the temperature schedule. The two- 
level system will be alternately exposed to a hot and a cold reservoir, which means that 
the function /?(t) in equation ([3]) will be a piecewise constant periodic function. During 
the first branch, it assumes the value during the second branch it attains the value 

This completes the description of the model. Any quantity describing the engine's 
performance can only depend on the parameters hi, /12, t±, and u. In the following 
we will focus on the characterization of the limit cycle, which the engine will approach 
at long times after a transient period. 
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(12) 



We start from the general solution (jl]) of the master equation ([T]). Owing to the 
Chapman- Kolmogorov condition ([6]), the propagator within the cycle is 

^_(t)M+(t+), t G [t+,tp] . 

Here the matrixes M±(t) evolve the state vector within the respective branches and have 
the form 

1 / 1 -1 \ . ... 1 / -1 -1 



,1 '\ [1 - exp (-z/t)] + i 



1 1 



where 



2 



dr exp \—v{t — r)] tanh < /3. 



dr exp \—v{t — r)] tanh < /5_ 



hi H ^ r 



(13) 



(14) 



(15) 



(16) 



Notice, both propagators M+(t) and ]R_(t) are given by the generic propagator (jlj). In 
order to get ffi+(t), we replace in equation ([5]) the initial position of the first energy h 
by hi, the driving velocity v by v+ = {h2 — hi)/t^, and we set t' = 0. Analogously, 
the propagator M_(t) follows from the generic propagator, if we replace h by h2, v by 
V- = {hi — h2)/t-, and t' by ^4.. 

The system state probabilities at the ends of the periods form a Markov chain and 
we are interested in its fixed point behavior. If we take the stationary state as the initial 
condition, the system revisits this special state at the end of the limit cycle. Therefore 
it suffices to solve the eigenvalue problem ]R_(t_)R_|_(t_|_) = |pstat^_ Solving the 

algebraic equation, the fixed point probabilities pf at the beginning (or end) of the 
limit cycle are 



Stat 
Pi 



1 - 



stat 



1 - 



e+(t+)exp(-z/t-)+e-(t- 
1 — exp {—i^tp) 



[17) 



These probabilities, and hence also the specific form of the limit cycle, depend solely on 
the model parameters. 

We now put aside the transitory regime and we focus entirely on the limit cycle. 
Generally speaking, the parametric plot of the occupation difference p{t) = pi{t) —p2{t) 
(the response) versus the energy of the first level Ei{t) (the driving) exhibits two 
possible forms which are exemplified in figure 1. First, we have a one- loop form 
which is oriented either clockwise or anticlockwise. For clockwise orientation, the work 
done by the engine on the environment during the limit cycle is negative, while for 
counter-clockwise orientation it is positive. Secondly, we can obtain a two-loops shape 
exhibiting again either positive or negative work on the environment. Slowing down 
the driving, the branches gradually approach the corresponding equilibrium isotherms 
p±{E) = — tanh(/5±_E/2). We postpone the further discussion to the section HI 
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Figure 1. The limit cycle for the two-stroke engine. The three graphs in the upper 
panel illustrate the case where /12 > > and the energy levels do not cross during 
their driving. On the left side we show E{t) = Ei{t) = —E2{t) and the response 
p{t) — pi{t) ~ p2{t). On the right hand side the parametric plot of the limit cycle 
in the p—E plane is displayed. The cycle starts in the upper vertex and proceeds 
counterclockwise, c.f. the arrows. The dashed and the dot-dashed curves show the 
equilibrium isotherms corresponding to the baths during the first and the second 
stroke, respectively. The parameters are: hi — 1 3, h2 = 5 J, t+ = 5 s, = 15 s, 
/3+ — 0.5 J~^, /3- = 0.1 J~^, V = 1 . The three graphs in the lower panel depict 
the case where /ii < < and the energies cross twice during the cycle. Except 
hi — —2.5 J, all parameters are as above. 

3. Probability densities for work and heat 

Heuristically, the underlying time-inhomogeneous Markov process D(t) can be conceived 
as an ensemble of individual realizations (sample paths). A realization is specified by 
a succession of transitions between the two states. If we know the number n of the 
transitions during a path and the times {tk}^=i at which they occur, we can calculate 
the probability that this specific path will be generated. A given paths yields a unique 
value of the microscopic work done on the system. For example, if the system is known 
to remain during the time interval [tfc, t/c-(-i], tfc+i > tk-, in the ith state, the work done 
on the system during this time interval is simply Ei{tk+i) — Ei(tk). Accordingly the 
probability of the paths gives the probability of the work. Viewed in this way, the 
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work itself is a stochastic process and we denote it as \N(t). We are interested in its 
probability density p{w, t) = { 6{\N{t) — w)), where (...) denotes an average over all 
possible paths. 

As a technical tool we introduce the augmented process {W(t), D(t)} which 
simultaneously reflects both the work variable and the state variable of the Markov 
process. The augmented process is again a time non-homogeneous Markov process. 
Actually, if we know at a fixed time t' both the present state variable j and the 
work variable w', then the subsequent probabilistic evolution of the state and work 
is completely determined. The work done during a time period [t',t], where t > t', 
simply adds to the present work w' and it only depends on the succession of the states 
after the time t'. And this succession by itself, as we know from section [21 cannot depend 
on the dynamics before time t'. 

The one-time properties of the augmented process will be described by the functions 
Prob { W(t) e{w,w + e) and D(t) = i \ W(t') = w' and D(t') = j } 



Gij{w, 1 1 w' ^ t') = lira 



(18) 

where i,j = 1,2. We represent them as the matrix elements of a single two-by-two 
matrix G{w, w'; t, t'), 

Gij{w, t\w',t') = {i\ G{w, 1 1 w\ t') I J ) . (19) 

Using this matrix notation, the Chapman- Kolmogorov condition for the augmented 
process assumes the form 

G(w, 1 1 w\ t') = I dw" G{w, 1 1 w", t") G(w", t" I w\ t') . (20) 



Here the matrix multiplication on the right hand side amounts for the summation 
over the intermediate states at the time t", and the integration runs over all possible 
intermediate values of the work variable w". The equation is valid for any intermediate 
time t" G [t',t]. Similarly to the preceding section [2l the Chapman-Kolmogorov 
condition can be used to connect two different propagators describing the time evolution 
of the augmented process within two branches of the driving cycle. Before we address 
this point, we focus on the generic situation. 

We need an equation which controls the time dependence of the propagator 
G{w,t\w',t') and which plays the same role as the rate equation ([1]) in the case of 
the simple two-state process. It reads 

(21) 

where the initial condition is G{w, t' \ w', t') = 6{w — w')I. This is a hyperbolic system 
of four coupled partial differential equations with time-dependent coefficients. It can 
be derived in several ways. For example, as explained in reference [2^, one considers 
at the time t the family of all realizations, which display at that time the work in 
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the infinitesimal interval {w, w + dw) and, simultaneously, which occupy a given state. 
During the infinitesimal time interval {t, t+dt), the number of such paths can change due 
to two reasons. First, while residing in the given state, some paths enter (leave) the set, 
because the energy levels move and an additional work has been done. Secondly, some 
paths can enter (leave) the described family because they jump out of (into) the specified 
state. These two contributions correspond to the two terms on the right hand side of 
equation fl2T|) . Another derivation [21j is based on an explicit probabilistic construction 
of all possible paths and their respective probabilities. 

Similar reasoning holds for the random variable Q(t) describing the heat accepted 
by the system from the environment, and for the internal energy U(t). The variable 
Q(t) is described by the propagator K(g,t | q',t') with the matrix elements 

( ,| , y Prob {Q(t) G (g, g + 6) A D(t) = i \ Q{f) = g' A D{t') = j} 
Kij{q,t \ q ,t) =\im . (22) 

e^O e 

It turns out that there exists a simple connection between the heat propagator and the 
work propagator G(w, t \ w', f). Since for each path, heat q and work w are connected 
by the first law of thermodynamics, we have q = Ei{t) — Ej{t') — w for any path which 
has started at time t' in the state i and which has been found at time t in the state j. 
Accordingly, 

nq,t\q',t')=( 9niun{t,t') -q,t\q' t') g,,Mt,t') -q,t\q',^ \ 

\g2i{u2i{t,t')-q,t\q',t') g22{u22{t,t')^q,t\q',t') j ' 

where Uij(t, t') = Ei{t)—Ej{t'). This relation can be written in the form of the symmetry 
relation 

G,j{uij{t, t')/2 + g, 1 1 g', t') = K,j{ui,{t, t')/2 -q,t\ q', t') . (24) 
3.1. Generic case-linear driving protocol 

For the linear driving protocol Ei[t) = h + v{t — t') = —E2{t) the first term in the curly 
brackets in equation ( 12T|) is time- independent. As for the second term, we use again the 
Glauber rates ([7]). Thereby the evolution equation (pT]) assumes the form 



(25) 

u f 1 -cexp[-n{t -t')]\] , , , 

+ l + cexp[-Qit-t')] [ -1 cexp[-n{t-t')] JjM^'^l^'^)' 

with the parameters c and Q introduced in connection with equation ([7j). 

We shall now employ the method described in reference by taking the double 
Laplace transformation with respect to the variables t and w. As shown in reference |21j . 
a special difference equation results, which can be solved exactly. Moreover, it is possible 
to carry out the final double inverse Laplace transformation. However, in [21], only the 
case Ei{0) = £^2(0) = has been studied. In the present context we need the solution 
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for a general initial energy difference 2h. It turns out that such generalization represents 
a nontrivial task. The constant c cannot be simply scaled off because it enters only the 
second term on the right hand side of equation f l2Tl) . In order to overcome this difficulty, 
we had to modify the procedure from reference j^T]. However, in view of the specific 
topic of the present paper, we refrain from giving the technical details and we proceed 
with the description of the final result. 

For the presentation of the result it is convenient to introduce the reduced work 
variable r] = rj{w, w') = 2[3{w — w') and the reduced time variable r = r(t, t') = VL{t — t'). 
Moreover, it is helpful to use the abbreviations 

l-x 1 



X 



exp 



7] 



y = exp 



— c ■ 



y_ 

1 + cx I + cy 



(26) 



For t> > 0, the result is 



[1 + c) exp(-r) 



ac 



(r - r/) + 0(r + r/)e(r - 7]) - x)y 



X 



1 + cexp(— r) 

2Fi(l + a, -a; 1; 0) , \n i ^ 2Fi(2 + a, 1 - a; 2; 0) 
-— — hl + al + cl + cxy)- — — z — 



1 



n f \ ' '\ ^r^^ ^ \r^^ \ « 2Fi(a, 1 - a; 1; 0) 
-GMV, - h , r ) = -0(r + ,)e(r - ,) acx y ^ ^^y^^ ^ ^^y.. 



1 



1 



— G2i(77, r I r') = -Q{r + r^)&{r - r/) ax'^ 



2Fi(l + a, -a; 1; 0) 
1 + cx)i+'*(l + cy)-'' 



(27) 
(28) 
(29) 



^^22(^7, r I r7',r') 
2Fi(a, 1 



1 + cexp(— r) 



1 + c 



ac 



{T + r]) + e(r + r])e{T - r]) —x''{l - y) 



X 



+ 



a; 1;. 



- (1 - a)(l +c)(l + cxy) 



2Fi(l + a, 2 - a; 2; , 



(30) 



[1 + cx)^+''{l + cy)^-"" (1 + cx)2+'^(l + c?/)2-«_ 

Here 5{-) is the Dirac delta-function, and B(-) is the Heaviside unit step function. The 
solution for f < follows from interchanging the indices 1 and 2 in equations fl271) - 
( !30|) . If /i = 0, then c = 1, and our results coincide with the formulae (49)-(52) in 
reference |21j . 



3.2. Piecewise linear periodic driving 

The generic result (!271) - (!30l) immediately yields the work and heat propagators for the 
individual branches in the protocol according to equation ffTTl) . We simply carry out the 
replacements described in the text following equation f|T5|) . We denote the corresponding 
matrices as G±{w,w' ,t) and K±{w,w',t). Then the Chapman- Kolmogorov condition 
( !20|) yields the propagator 

G+{w,o,t), te[o,t+], 



ip(:w,t) 



h2—hi 



-{h2-hi) 



dw'G^{w,w',t)G+{w',0,t+) , te [t+,tp] 



(31) 
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As demonstrated above, the heat propagator Kp{w,t) for the hmit cycle is connected 
with the work propagator Gp{w,t) through simple shifts of the independent variable 
w. Specifically, we get ( i \ Kp(g, t)\j) = {i\ Gp{uij{t) — q,t)\j) with U2i{t) = —uuit), 
U22{t) = -■ uii(t), and 

Unit) = ^^^^t, m2{t) = 2h + ^^^^t, t G [0,t+], (32) 

Unit) = {h2 - h) (^1 - ^j^) , ui2{t) = h2 + h- ^'~^' (t - U), t e [U, tp]. (33) 

In the last step we take into account the initial condition (fTTl) at the beginning of 
the limit cycle and we sum over the final states of the process D{t). Then the probability 
density for the work done on the system during the limit cycle reads 

2 

Pp{w,t) = J2{^\Gpiw,t)\f'^') . (34) 

i=l 

Similarly, the probability density for the head accepted within the limit cycle is 

2 

xA1,t) = Y,{^\Kp{q,t)\f'^') . (35) 

i=l 

These two functions represent the main results of the present Section. They are 
illustrated in figures 2-4. We discuss their main features in section O 



4. Engine performance 

As shown in section [2l the occupation probabilities during the limit cycle are Mp(t) 
with Mp(t) given by equation ( JT2ll . These probabilities are ensemble averaged quantities 
and cannot describe fluctuations of the engine's performance. But they render the 
energetics in terms of mean values as we discuss now. 

During the limit cycle, the internal energy U(t) = Yl'i=i Ei(t)pi(t) changes as 

2 2 

fu{t) = J2E^{t)f^P^it) + J2P^it)^t^^it) = ^ [QW + W^W] , t ^ [0, tp] . (36) 

i=l i=l 

Here Q{t) = (Q(t) ) is the mean heat received from the reservoirs during the period 
between the beginning of the limit cycle and the time t. Analogously W{t) = ( W(t) ) is 
the mean work done on the system from the beginning of the limit cycle till the time t. If 
W{t) < 0, the positive work —W{t) is done by the system on the environment. Therefore 
the oriented areas enclosed by the limit cycle in figure 1 and in figure 4 represent the 
work Wout = ~W{tp) done by the engine on the environment per cycle. These areas 
approach maximal absolute values in the quasi-static limit. The internal energy, being 
a state function, fulfills U(tp) = U{0). Therefore, if the work PKut is positive, the same 
total amount of heat has been transferred from the two reservoirs during the cycle. The 
case Wout > cannot occur if both reservoirs would have the same temperature. That 
the perpetuum mobile is actually forbidden can be traced back to the detailed balance 
condition in ([T]). 
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Figure 2. The probability density pp{w,t) as a function of the work w for the same 
parameters as in figure 1 (with positive hi): a) t = ^t^ (middle of the first stroke), b) 
t = t+ (end of the first stroke), c) t = t+ + ^t^ (middle of the second stroke), and d) 
t = t+ + t- (end of the limit cycle). The triangle on the work axis marks the mean 
work W{t) at the corresponding times. The singular parts of pp{w,t) are marked by 
arrows, where the arrow heights equal the weights of the corresponding delta functions 
[for example, in panel a), the left arrow height gives the probability that the system is 
initially in the second state and remains in it between the beginning of the cycle and 
the time t = ^t+; then the work done on the system equals — |(/i2 — hi)]. 



We denote the system entropy at time t as Ss{t), and the reservoir entropy at time 
t as Sj-it). They are given by 

^'^^^ --[p,it)\np^{t)+,it)\np,it)], (37) 



B 



k 

S (t) /"*+ d /"^p d 

-J^ = -(3^l dt'Ei{t')-[p,{t')-p,{t')]-(3.j^ dt'E,{t')-[pi{t')-p,{t')]. (38) 

Upon completing the cycle, the system entropy re-assumes its value at the beginning of 
the cycle. On the other hand, the reservoir entropy is controlled by the heat exchange. 
Owing to the inherent irreversibility of the cycle we observe always a positive entropy 
production per cycle, Sj.(tp) — 5*1.(0) > 0. The total entropy 5'tot(^) = Ss(t) + S^(t) 
increases for any t G [0,tp]. The rate of the increase is the larger the stronger is the 
representative point in the p—E diagram deviates from the corresponding equilibrium 
isotherm (a strong deviation, e.g., can be seen in the p—E diagram in figure 4c). Due to 
the instantaneous exchange of the baths at times t+ and t^ + t_ in the model considered 
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Figure 3. The probability density Xp{l: t) a function of the heat q and for the same 
parameters as in figme 1 (with positive hi): a) t = ^i-)- (middle of the first stroke), 
h) t = (end of the first stroke), c) t = t+ + it_ (middle of the second stroke), and 
d) < = + t_ (end of the limit cycle). The triangles on the heat axis mark the mean 
heat Q{t) at the corresponding times. The singular parts of Xp('Z: t) are marked by the 
arrow, where the arrow height equals the weight of the corresponding delta function. 
For example, in a), the height of the arrow gives the probability that there was no 
transition between the states from the beginning of the cycle till the observation time 
t = The heat exchanged in this case is zero. 



here, a strong increase of Stot{t) always occurs after these time instants. A representative 
example of the overall behavior of the thermodynamic quantities (mean work and heat, 
and entropies) during the limit cycle is shown in figure 5. 

An important characteristics of the engine is its power output Pout and its efficiency 
/i. They are defined as 

^out = — — , /i = , (39) 

t-p ^5 in 

where Qi^ is the total heat absorbed by the system per cycle. The performance of the 
engine characterized by the output work, efficiency, output power, and entropies from 
equations (^71) and are shown in figure 6 and figure 7. 

In figure 6 the performance is displayed as a function of the cycle duration tp for 
t+ = t_ = tp/2. With increasing tp, the output work and the efficiency increase whereas 
the output power and the entropy production first increase up to a maximum and 
thereafter they decrease when approaching the quasi-static limit (tp — > oo). Notice that 
the maximum efficiency and output power occur at different values of tp. In figure 6a) 
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Figure 4. Probability densities pp{w,tp) and Xp(9j^p) for the work and heat for 
four representative sets of the engine parameters. For every set we show also the limit 
cycle in the p—E plane, where the corresponding equilibrium isotherms are marked 
by dashed (first stroke) and dot-dashed (second stroke) lines. In all cases we choose 
hi ~ 1 3, h2 — 5 3, and 1^ = 1 s~^. The remaining parameters are a) t+ = 50 s, 
<_ = 10 s, /3+ = 0.5 J^^, /9_ = 0.1 J^^ (bath of the first stroke is colder than of the 
second stroke), b) t+ = 50 s, t- = 10 s, (3+ = 0.1 J~^, /3_ = 0.5 (exchange of /3+ 
and P- as compared to case a), leading to a change of the traversing of the cycle from 
counter-clockwise to clockwise and a sign reversal of the mean values W{tp) = { \N{tp) ) 
and Q{tp) = {Q{tp) )), c) t+ = 2 s, t_ = 2 s, /3+ = 0.2 J^^, /3_ = 0.1 J^^ (a strongly 
irreversible cycle traversed clockwise with positive work), d) t+ — 20 s, t- = 1 s, 
f3± = 0.1 J^^ (no change in temperatures, but large difference in duration of the two 
strokes; W{tp) is necessarily positive). 



we show also the standard deviation of the output work, which was calculated from the 
work probability density Pp{w,tp). Finally, let us note that the values /?+ = 0.5 J^^ 
and /3_ = 0.1 used in figure 6 give the Carnot efficiency /ic = 0.8. This should be 
compared with the efficiency of the engine for a long period tp, that is, with the value 
fi ^ 0.6. As discussed above, the Carnot efficiency cannot be reached here even for 
tp oo, due to the immediate temperature changes at times t+ and t+ + 1_. 

In figure 7 we have fixed tp and plotted the behavior as function of the time 
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Figure 5. Thermodynamic quantities as functions of time during the hmit cycle 
for the same set of parameters as in the upper panel of figure 1 (positive hi). Left 
panel: internal energy, mean work done on the system, and mean heat received from 
both reservoirs; the final position of the mean work curve marks the work done on 
the system per cycle W{tp). Since W{tp) < 0, the work Wout — —W{tp) has been 
done on the environment. The internal energy returns to its original value and, after 
completion of the cycle, the absorbed heat Q{tp) equals the negative work —W{tp). 
Right panel: entropy Ss{t) of the system and Sr{t) of the bath, and their sum Stot{t); 
after completing the cycle, the system entropy re-assumes its initial value. The 
difference S'tot(ip) — 5*101(0) > equals the entropy production per cycle. It is always 
positive and quantifies the degree of irreversibility of the cycle. 

asymmetry (or time splitting) parameter A = (t+ — tJ)/tp. As can be seen from 
the upper three panels in figure 7, there exist also a maximal efficiency and a maximal 
output power with respect to a variation of the time asymmetry parameter (as long as 
the engine performs work, i.e., Wont > 0). Again, the optimal parameter A, where these 
maxima occur, is different for the efficiency and output power. In a reversed situation, 
considered in the lower three panels in figure 7, where the work is performed on the 
engine (Wout < 0), minima of the efficiency and output power occur. 

5. Discussion 

The overall properties of the engine critically depend on the two dimensionless 
parameters a± = i^/{'2f3±\v±\). We call them reversibility parameters. For a given 
branch, say the first one, the parameter a+ represents the ratio of two characteristic 
time scales. The first one, l/iy, is given by the attempt rate of the internal transitions 
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Figure 6. The engine performance versus the duration of the limit cycle t± for 
= = itp and otherwise the same parameters as in the upper part of figure 1 
(positive hi). Both the output work Wout in a) and the efficiency fi in b) increase with 
tp. The output power Pout in c) assumes a maximum at a special cycle duration. The 
dashed line in a) marks the standard deviation of the output work, calculated from the 
work density pp(w, ip). Notice that the work fiuctuation is comparatively high close to 
the cycle duration, where the maximal output power is found. In the long-period limit 
tp — > cxD, the cycle still represents a non-equilibrium process (due to the construction of 
the model, see text), and hence the entropy production S'tot(tp) — 5*101(0) in d) remains 
positive, approaching a specific asymptotic value. 

[29] . The second scale is proportional to the reciprocal driving velocity. Contrary to the 
first scale, the second one is fully under the external control. Moreover, the reversibility 
parameter is proportional to the absolute temperature of the heat bath. 

Let us first consider the work probability density (IMIl within the first stroke in the 
case h2 > hi, c.f. figure. 2a). In essence, pp{w,t) is given by a linear combination of 
the functions equations fl27|) - fl30|) . It vanishes outside the common support [—v+t,v+t] 
which broadens linearly in time. Besides the continuous part located within the support, 
the diagonal elements ga^w, 0, t, 0) display a singular part represented by delta functions 
at the borders of the support. The delta functions correspond to the paths with no 
transitions between the states. Specifically, the weight of the delta function located at 
w = v^t represents the probability that the system starts in the first state and remains 
there up to time t. The weight corresponding to the first level decreases with increasing 
time and vanishes for t — > oo. On the contrary, the weight of the delta function at — v+i 
approaches the nonzero limit 2/3_|_/(l -f c+)°+ for t oo, which is the probability that 
the a path starts in the second state and never leaves it. 

Within the second stroke, the density Pp{w,t) results from the integral of the 
propagators for the individual strokes, c.f. equation ( l3T|) . Due to the integration, the 
singular parts of the cycle propagator G+(w, 0; t) are now situated inside the support, at 
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Figure 7. The engine performance characterized by the efficiency //, output power 
Pout, and entropy production Stot{tp) — Stot{0), as a function of the asymmetry 
parameter A = — t^)/tp for a fixed period ip = 20 s and the same parameters 
hi = 1 3, h2 — 5 3, v — 1 s~^ as in the upper panel of figure 1. In al)-cl) the 
bath during the first stroke is colder than during the second stroke: (3+ = 0.5 and 
/3_ = 0.1 J-^ Notice that the value A of maximum efficiency does not correspond 
to that of maximum output power. In a2)-c2) the reciprocal bath temperatures are 
interchanged compared to cases al)-cl), (3+ = 0.1 and /3_ = 0.5 J~^. The dashed 
curves in bl) and b2) show the standard deviation of the output power calculated from 

Pp{w,tp). 



the values w = —v+t+ + \v^\(t — t+) and w = v^t+ — \v^\{t — t+). The two delta functions 
approach each other and, upon completing the cycle, they coincide at the point w = 0. 
The nonsingular component of the density is no more continuous |30]. The jumps are 
located at the positions of the delta functions and their magnitudes correspond to the 
weights of the delta functions (for a discussion of the origin of these jumps, see [3T]). 

If both reversibility parameters a± are small, the isothermal processes during both 
branches strongly differs from the equilibrium ones. The signature of this case is a 
flat continuous component of the density pp{w,t) and a well pronounced singular part. 
The strongly irreversible dynamics occurs if one or more of the following conditions 
hold. First, if u is small, the transitions are rare and the occupation probabilities 
of the individual energy levels are effectively frozen during long periods of time. 
Therefore they lag behind the Boltzmann distribution which would correspond to the 
instantaneous positions of the energy levels. More precisely, the population of the 
ascending (descending) energy level is larger (smaller) than it would be during the 
corresponding reversible process. As a result, the mean work done on the system is 
necessarily larger than the equilibrium work. Secondly, a similar situation occurs for 
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large driving velocities v±. Due to the rapid motion of the energy levels, the occupation 
probabilities again lag behind the equilibrium ones. Thirdly, the strong irreversibility 
occurs also in the low temperature limit. In the limit a± 0, the continuous part 
vanishes and Pp{w,tp) = 6{w). 

In the opposite case of large reversibility parameters a±, both branches in the 
p—E plane are located close to the reversible isotherms. The singular part of the 
density Pp{w,t) is suppressed and the continuous part exhibits a well pronounced peak. 
From general considerations [H], the density must approach a Gaussian shape. Our 
results allow a detailed study of this approach. Let us denote as F{(3, E) the free 
energy of a two level system with energies ±-E at temperature T = 1/(^3/9, i.e., 
F{I3,E) = -iln[2cosh(/3E)]. Let us further define 



F{P_,E,{t)) - F{P_,E,{U)) + F{P+,E,{U)) - F{P+,E,{0)), t G [t+,tp]. 

(40) 



This is simply the reversible work done on the system if we transform its state from 
the initial equilibrium state (with the energies fixed at ±£"1(0)) to another equilibrium 
state (with the energies fixed at the values ±Ei{t)). For large reversibility parameters 
a±, the peak of the work density Pp{w, tp) occurs in the vicinity of the value Wrevit) and 
with increasing a±, the peak collapses to a delta function. 



The main features of the heat probability density Xpil^t) from equation (1351) are, 
as we have seen in section [3], closely related to the work through simple shifts of the 
independent variable q. However, there are some interesting differences. While the work 
is conditioned by the external driving, the heat exchange occurs as a consequence of the 
transitions between the system states. The instantaneous positions of the energies at the 
instant of the transition give the magnitude of the heat exchange related with the given 
transition. From this perspective, if there are no transitions, the exchanged heat is zero. 
As a consequence, the singular part of the probability density Xpily t) is always situated 
at g = and the weight of the delta function at origin equals the sum of the weights 
of the delta functions in the work density pp{w,t). The support of the heat density is 
given by the largest possible value of the level splitting during the limit cycle. Within 
the first stroke the support broadens linearly with time as [—2hi — 2v+t, 2hi + 2v+t], up 
to its maximum width [—2/12,2/12] at the end of the stroke. Within the second stroke 
the energy difference decreases and the support remains unchanged. The non-singular 
part of the heat density always displays discontinuities inside the support, even during 
the first stroke. In contrast to ( z | Gp{w, t)\j), the individual elements ( i \ Kp(g, t)\j) 
in equation fl35l) have different supports. 

In the the strongly reversible regime each element {i\Gp{w,t)\j) exhibits a 
Gaussian shape situated at W^ev{t)- The shift transformation maps the Gaussian 




tG [0,t+], 



lim Pp{w,t) = 6{W - Wrevit)) . 



(41) 
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function onto four different positions depending on the specific matrix element 
{i I Kp(g, t)\j) in question. In the reversible limit we have 



Using this form in equation ( l35i) and calculating the mean accepted heat, we get 
Q{t) = U(t) — U{0) — Wj-ev(t). In the opposite limit, if a± — > 0, we have Xpil^t) ~^ ^il) 
for any t. 

According to the second law of thermodynamics, the mean work W{t) = ( W(t) ) 
must fulfill > |VFrcv(i)|- On the other hand, there always exists a fraction of 

the paths which, individually, display the inequality |w(path, t)| < |lVrev(^)|, where 
w(path, t) denotes the work done on the system if it evolves along the indicated path. 
Using the exact work probability density, we can calculate the total weight of these 
trajectories. Specifically, in the case W^evit) > 0, 



If Wrevit) < 0, we would have to integrate over the interval (Wrevit), co, )• 

Let us finally note that in view of the rather complex structure of the work and 
heat probability densities, we performed several independent tests. First of all, the 
densities Pp{w,t) and XpiQ^t) must be nonnegative functions fulfilling the normalization 
conditions, e.g. J^^dwpp{w,t) = 1 for any t G [0,tp]. Secondly, we have two 
different procedures to calculate the first moment W{t) = (W(t)). One can either 
start with the density pp{w,t) and evaluate the required w-integral, or one directly 
employs the solution of the rate equation as in Sec. IV. Another inspection is based 
on the Jarzynski identity [11 [8]. In our setting, consider the case /3± = /3. After 
completing the cycle, the system returns to the original state. Therefore we have 
W^rev(^p) = F{P,Ei{tp)) — F{P, Ei{0)) = and the Jarzynski identity reduces to 
(exp [—(3\N{tp)] ) = 1. Using the explicit form of the work probability density we have 
verified that the integral dw exp{—Pw)pp{w,tp) actually equals one. Finally, we 
have studied the probability densities pp{w,t), Xpilyt) by computer simulation. In 
fact, we have developed two exact simulation methods. Each of them uses a specific 
algorithm to generate paths of the time-non-homogeneous Markov process D(t). Parts 
of these simulation results have been published in [STj and confirm the analytical results. 

6. Conclusions 

We have investigated a simple example of a microscopic heat engine, which is 
exactly solvable. Based on mean thermodynamic quantities, the engine performance 
is characterized by the occupation probabilities of the energy levels following from 
the master equation. The more challenging exact calculation of the work and heat 
probability densities allowed us to study the fluctuation properties in detail. A 
notable result is that the engine can be tuned to maximize its output power, but the 



lim ( z I Kp{q, t)\j) = 6{q- Uij{t) + Wre^t) ) 



(42) 



a± 




(43) 
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fluctuations of tliis quantity in tlie corresponding optimal regime of control parameters 
are comparatively high. 

The present setting can be expanded in various directions. One can address various 
problems concerning the thermodynamic optimization. Another option would be the 
embodiment of additional (e.g., adiabatic) branches. The role of the working medium 
can be assigned to other systems that exhibit more complicated dynamics (e.g., diffusing 
particles in the presence of time-dependent forces, or, variants of the generalized master 
equation). It would be also interesting to investigate settings with a nonlinear driving 
of the energy levels. A nontrivial generalization would be the inclusion of a third energy 
level. Having the three levels one can couple the system (different pairs of forth-back 
transitions between the levels) simultaneously to reservoirs at different temperatures, 
so that the system approaches a non-equilibrium steady state without driving |32j . 
Including a driving and forming an operational cycle, there is no serious obstacle in 
repeating the present analysis for this system, which has some additional intriguing 
properties compared to the two-level system considered here (as, e.g., negative specific 
heats) . 

Another possibility is an incorporation of specific forms of transition rates [33] that 
describe the stretching of biomolecules in some realistic manner. In such problem, the 
histogram of the work is experimentally accessible [33]. Particularly, in the experiments 
one can also determine the probability of having certain number of transitions between 
the folded and the unfolded conformation of the biomolecule during its mechanical 
stretching [33]. In our formulation, this information is encoded in the counting statistics 
of the underlying random point process [3^ and can be extracted from the perturbation 
expansion of the propagators which solve our dynamical equations. Calculations in this 
direction are in progress and will be reported elsewhere. 
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